I'm having a problem with a few different version of my run_mcmc function, where it's giving shitty contours. I'm gonna put it into a notebook so I can more easily make plots and step through different processes in it.
In [1]:
from pearce.emulator import OriginalRecipe, ExtraCrispy
from pearce.mocks.customHODModels import *
from pearce.mocks import cat_dict
import numpy as np
from os import path
In [2]:
from multiprocessing import cpu_count
import warnings
from itertools import izip
import numpy as np
import emcee as mc
from scipy.linalg import inv
In [3]:
from matplotlib import pyplot as plt
%matplotlib inline
import seaborn as sns
sns.set()
In [4]:
def lnprior(theta, param_names, *args):
for p, t in izip(param_names, theta):
low, high = _emu.get_param_bounds(p)
if np.isnan(t) or t < low or t > high:
return -np.inf
return 0
In [5]:
def lnlike(theta, param_names, fixed_params, r_bin_centers, y, combined_inv_cov, obs_nd, obs_nd_err, nd_func_name):
param_dict = dict(izip(param_names, theta))
param_dict.update(fixed_params)
y_bar = _emu.emulate_wrt_r(param_dict, r_bin_centers)[0]
# should chi2 be calculated in log or linear?
# answer: the user is responsible for taking the log before it comes here.
delta = y_bar - y
print 'y',y
print 'ybar',y_bar
#print y_bar
chi2 = -0.5 * np.dot(delta, np.dot(combined_inv_cov, delta))
return chi2# - 0.5 * ((obs_nd - getattr(_cat, nd_func_name)(param_dict)) / obs_nd_err) ** 2
In [6]:
def lnprob(theta, *args):
lp = lnprior(theta, *args)
if not np.isfinite(lp):
return -np.inf
return lp + lnlike(theta, *args)
In [7]:
def _run_tests(y, cov, r_bin_centers, param_names, fixed_params, ncores, nd_func_name):
assert ncores == 'all' or ncores > 0
if type(ncores) is not str:
assert int(ncores) == ncores
max_cores = cpu_count()
if ncores == 'all':
ncores = max_cores
elif ncores > max_cores:
warnings.warn('ncores invalid. Changing from %d to maximum %d.' % (ncores, max_cores))
ncores = max_cores
# else, we're good!
assert y.shape[0] == cov.shape[0] and cov.shape[1] == cov.shape[0]
assert y.shape[0] == r_bin_centers.shape[0]
# check we've defined all necessary params
assert _emu.emulator_ndim == len(fixed_params) + len(param_names) + 1 # for r
tmp = param_names[:]
assert not any([key in param_names for key in fixed_params]) # param names can't include the
tmp.extend(fixed_params.keys())
assert _emu.check_param_names(tmp, ignore=['r'])
assert hasattr(_cat, nd_func_name)
return ncores
In [8]:
def _resume_from_previous(resume_from_previous, nwalkers, num_params):
# load a previous chain
# TODO add error messages here
old_chain = np.loadtxt(resume_from_previous)
if len(old_chain.shape) == 2:
c = old_chain.reshape((nwalkers, -1, num_params))
pos0 = c[:, -1, :]
else: # 3
pos0 = old_chain[:, -1, :]
return pos0
In [9]:
def _random_initial_guess(param_names, nwalkers, num_params):
"""
Create a random initial guess for the sampler. Creates a 3-sigma gaussian ball around the center of the prior space.
:param param_names:
The names of the parameters in the emulator
:param nwalkers:
Number of walkers to initiate. Must be the same as in resume_from_previous
:param num_params:
Number of params to initiate, must be the same as in resume_from_previous
:return: pos0, the initial position of each walker for the chain.
"""
pos0 = np.zeros((nwalkers, num_params))
for idx, pname in enumerate(param_names):
low, high = _emu.get_param_bounds(pname)
pos0[:, idx] = np.random.randn(nwalkers) * (np.abs(high - low) / 6.0) + (low + high) / 2.0
# TODO variable with of the initial guess
return pos0
In [10]:
def run_mcmc(emu, cat, param_names, y, cov, r_bin_centers, obs_nd, obs_nd_err, nd_func_name, \
fixed_params={}, resume_from_previous=None, nwalkers=1000, nsteps=100, nburn=20, ncores='all'):
_emu = emu
_cat = cat
global _emu
global _cat
ncores= _run_tests(y, cov, r_bin_centers,param_names, fixed_params, ncores, nd_func_name)
num_params = len(param_names)
combined_inv_cov = inv(_emu.ycov + cov)
sampler = mc.EnsembleSampler(nwalkers, num_params, lnprob,
threads=ncores, args=(param_names, fixed_params, r_bin_centers, y, combined_inv_cov, \
obs_nd, obs_nd_err, nd_func_name))
if resume_from_previous is not None:
try:
assert nburn == 0
except AssertionError:
raise AssertionError("Cannot resume from previous chain with nburn != 0. Please change! ")
# load a previous chain
pos0 = _resume_from_previous(resume_from_previous, nwalkers, num_params)
else:
pos0 = _random_initial_guess(param_names, nwalkers, num_params)
return pos0, (param_names, fixed_params, r_bin_centers, y, combined_inv_cov, \
obs_nd, obs_nd_err, nd_func_name)
sampler.run_mcmc(pos0, nsteps)
chain = sampler.chain[:, nburn:, :].reshape((-1, num_params))
In [11]:
training_dir = '/u/ki/swmclau2/des/PearceLHC_wp_z_sham_emulator/'
em_method = 'gp'
split_method = 'random'
load_fixed_params = {'z':0.0}
In [12]:
emu = ExtraCrispy(training_dir,10, 2, split_method, method=em_method, fixed_params=load_fixed_params)
In [13]:
#Remember if training data is an LHC can't load a fixed set, do that after
fixed_params = {'f_c':1.0}#,'logM1': 13.8 }# 'z':0.0}
cosmo_params = {'simname':'chinchilla', 'Lbox':400.0, 'scale_factors':[1.0]}
cat = cat_dict[cosmo_params['simname']](**cosmo_params)#construct the specified catalog!
In [14]:
#mbc = np.loadtxt('/nfs/slac/g/ki/ki18/des/swmclau2/AB_tests/mbc.npy')
#cen_hod = np.loadtxt('/nfs/slac/g/ki/ki18/des/swmclau2/AB_tests/cen_hod.npy')
#sat_hod = np.loadtxt('/nfs/slac/g/ki/ki18/des/swmclau2/AB_tests/sat_hod.npy')
#cat.load_model(1.0, HOD=(HSAssembiasTabulatedCens, HSAssembiasTabulatedSats),\
# hod_kwargs = {'prim_haloprop_vals': mbc,
# 'cen_hod_vals':cen_hod,
# 'sat_hod_vals':sat_hod})
#cat.load_catalog(1.0)
cat.load(1.0, HOD='redMagic')
In [15]:
emulation_point = [('f_c', 0.2), ('logM0', 12.0), ('sigma_logM', 0.366),
('alpha', 1.083),('logM1', 13.7), ('logMmin', 12.233)]
#emulation_point = [('mean_occupation_centrals_assembias_param1',0.6),\
# ('mean_occupation_satellites_assembias_param1',-0.7)]
em_params = dict(emulation_point)
em_params.update(fixed_params)
#del em_params['z']
In [16]:
#rp_bins = np.logspace(-1.1,1.6,18)
#rp_bins.pop(1)
#rp_bins = np.array(rp_bins)
rp_bins = np.loadtxt('/nfs/slac/g/ki/ki18/des/swmclau2/AB_tests/rp_bins.npy')
rpoints = (rp_bins[1:]+rp_bins[:-1])/2.0
In [17]:
wp_vals = []
nds = []
for i in xrange(2):
cat.populate(em_params)
wp_vals.append(cat.calc_wp(rp_bins, 40))
nds.append(cat.calc_number_density())
In [18]:
#y = np.mean(np.log10(np.array(wp_vals)),axis = 0 )
y = np.log10(np.loadtxt('/nfs/slac/g/ki/ki18/des/swmclau2/AB_tests/sham_vpeak_shuffled_wp.npy'))
# TODO need a way to get a measurement cov for the shams
cov = np.cov(np.log10(np.array(wp_vals).T))#/np.sqrt(50)
#obs_nd = np.mean(np.array(nds))
obs_nd = np.loadtxt('/nfs/slac/g/ki/ki18/des/swmclau2/AB_tests/sham_vpeak_shuffled_nd.npy')
obs_nd_err = np.std(np.array(nds))
In [19]:
param_names = [k for k in em_params.iterkeys() if k not in fixed_params]
nwalkers = 10
nsteps = 10
nburn = 0
In [20]:
pos0, args = run_mcmc(emu, cat, param_names, y, cov, rpoints,obs_nd, obs_nd_err,'calc_analytic_nd', fixed_params = fixed_params,\
nwalkers = nwalkers, nsteps = nsteps, nburn = nburn)#,\
#resume_from_previous = '/u/ki/swmclau2/des/PearceMCMC/100_walkers_1000_steps_chain_shuffled_sham_no_nd.npy')#, ncores = 1)
In [21]:
from corner import corner
In [22]:
for t in pos0:
print lnlike(t, *args)
print
In [23]:
print cov
In [24]:
np.diag(emu.ycov)
Out[24]:
In [ ]: